library(SDMTools)#load necessary libraries
work.dir = 'E:/SEABat/'; setwd(work.dir) #define & set the working directory

#setup some plot parameters
legend.pnts = cbind(c(93,91,91,93),c(2,2,7.5,7.5))

threshold=0.335
current = read.asc(paste('fake.current.asc',sep=''))
current[which(is.finite(current) & current>threshold)] = 1 #set all parts of the species distribuiton to 1
current[which(current<1)]=0

cs.current=ClassStat(current, latlon=TRUE)
cs.current$total.area = cs.current$total.area / 1000000000 #area, 1000s km^2
	  
threshold=0.390
future = read.asc(paste('fake.future.asc',sep=''))
future[which(is.finite(future) & future>threshold)] = 2 #set all parts of the species distribuiton to 1
future[which(future<1)]=0

cs.future=ClassStat(future, latlon=TRUE)
cs.future$total.area = cs.future$total.area / 1000000000 #area, 1000s km^2

overlap=current+future

cs.overlap=ClassStat(overlap, latlon=TRUE)
cs.overlap$total.area = cs.overlap$total.area / 1000000000 #area, 1000s km^2

out=matrix(NA,nr=1,ncol=3)
out=as.data.frame(out)
colnames(out)=c('current','future','overlap')
out$current=cs.current[2,3]
out$future=cs.future[2,3]
out$overlap=cs.overlap[4,3]

write.csv(out,'area.csv')

cols = c('gray86','red3','forestgreen','blue')
png('overlaptest.png', width=nrow(bat.asc),height=ncol(bat.asc),, units='px', res=300, pointsize=5, bg='white') #start the plot
	par(mar=c(0,0,0,0)+0.1) #remove any plot margins
	image(overlap, ann=FALSE,axes=FALSE,col=cols, zlim=c(0,3)) #plot the richness
	legend.gradient(legend.pnts,cols=cols,limits=c(0,3), title='overlap', cex=2)
dev.off()
